%function hu_calibration(folder,dim,scans,ii,jj,kk)


function hu_fit=hu_calibration_fitting_toolbox(patient)
 
    %|------------------------------------------------------------------------|
    %| 
    %|  Then calculate linear fit each for slice to the deformed bellow
    %|  voltages;
    %------------------------------------------------------------------------
    %|      Dependancies;                                                     |
    %|                               |
    %------------------------------------------------------------------------
    %   This file is part of the
    %   5D-Novel4DCT Toolbox  ("Novel4DCT-Toolbox")
    %   DH Thomas, Ph.D
    %   University of California, Los Angeles
    %   Contact: mailto:dhthomas@mednet.ucla.edu
    %------------------------------------------------------------------------
    % $Author: DHThomas $	$Date: 2014/04/01 10:23:59 $	$Revision: 0.1 $


step=1;

ii=1:patient.dim(1);
jj=1:step:patient.dim(2);
kk=1:step:patient.dim(3);

[~,grid_z]=meshgrid(jj,kk);

grid_zs=permute(repmat(grid_z,[1,1,patient.scans]),[3,2,1]);
fpolyvalalt = @(p,x)sum(x*p,2)';

hu_fit=zeros(2,size(ii,2),size(jj,2),size(kk,2));

% if matlabpool('size')== 0 % checking to see if my pool is already open
%     matlabpool 4
% end

pctRunOnAll javaaddpath E:\ParforProgMonv2\java
progressStepSize = 1;
ppm = ParforProgMon('                      HU Correction ',size(ii,2),  progressStepSize,800, 100);

 tic
           pctRunOnAll  warning('off')

     parfor slice=ii;
         display(sprintf('Slice %d of %d',slice,size(ii,2)));
         
         [ dvf]=read_dvf_elastix_slice_toolbox(patient.dvf_folder,slice);
         
         img=load([patient.image_slice_folder sprintf('/deformed_images_slice_%d',slice)]);
         
         voxel_w=double(bsxfun(@plus,grid_zs,dvf.w(1:patient.scans,jj,kk)));
         voxel_w=voxel_w(1:patient.scans,:);
         img=img.img_tmp(1:patient.scans,jj,kk);
         img=img(1:patient.scans,:);
         hu_volt=zeros(patient.scans,size(voxel_w,2));
         
%          hu_volt_tmp=interp1(kk,bellows_volt_drifted(kk,:),voxel_w,'nearest');
%          hu_fit_tmp=zeros(2,size(jj,2)*size(kk,2));
%          for scan=1:scans;
%              hu_volt_tmp10=squeeze(hu_volt_tmp(scan,:,:))';
%              hu_volt(scan,:)=hu_volt_tmp10(scan,:);
%          end
         
         for scan=patient.ref:patient.scans;
             hu_volt(scan,:) = nakeinterp1(kk',patient.bellows_volt_drifted(kk,scan),voxel_w(scan,:)');            
         end
         
         
          hu_fit_tmp = zeros(2,size(img,2));
          tic
         for vox=1:size(voxel_w,2)
% tic
              hu_fit_tmp(:,vox)=[ones(patient.scans,1) hu_volt(:,vox)]\img(:,vox);
%              hu_fit_tmp_robust(:,vox) = robustfit((hu_volt(:,vox)),img(:,vox));
%              hu_fit_tmp_robust(:,vox) = regress((hu_volt(:,vox)),img(:,vox));
%              toc
         end
          toc
%          hu_fit_tmp1=reshape(hu_fit_tmp(1,:),size(jj,2),size(kk,2));
%          hu_fit_tmp2=reshape(hu_fit_tmp(1,:),size(jj,2),size(kk,2));
         
%         [hu_fit_smooth1 s1]=smoothn(hu_fit_tmp1,2.5, 'robust');
%         [hu_fit_smooth2 s2]=smoothn(hu_fit_tmp2,2.5,'robust');
% sigmaS=5;
% sigmaR=15;
% samS=5;
% samR=15;
% 
% %get the filter result
% Ibf=bilateral3(I, sigmaS,sigmaS,sigmaR,samS,samR);

         warning on
%          hu_volt=reshape(hu_volt,scans,size(jj,2),size(kk,2));
         hu_fit(:,slice,:,:)=reshape(hu_fit_tmp,2,size(jj,2),size(kk,2));
         
         ppm.increment();
     end
 toc

save([patient.folder sprintf('/hu_fit')],'hu_fit');
  


    